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ABSTRACT 

In order to investigate formation of relativistic jets at the center of a progenitor of a long gamma- 
ray burst (GRB), we develop a two-dimensional general relativistic magnetohydrodynamic (GRMHD) 
code. We show the code passes many, well-known test calculations, by which the reliability of the 
code is confirmed. Then we perform a numerical simulation of a collapsar using a realistic progenitor 
model. It is shown that a jet is launched from the center of the progenitor. We also find that the 
mass accretion rate after the launch of the jet shows rapid time variability that resembles to a typical 
time profile of a GRB. The structure of the jet is similar to the previous study: a poynting flux jet 
is surrounded by a funnel-wall jet. Even at the final stage of the simulation, bulk Lorentz factor of 
the jet is still low, and total energy of the jet is still as small as 10 48 erg. However, we find that the 
energy flux per unit rest-mass flux is as high as 10 2 at the bottom of the jet. Thus we conclude that 
the bulk Lorentz factor of the jet can be potentially high when it propagates outward. It is shown 
that the outgoing poynting flux exists at the horizon around the polar region, which proves that the 
Blandford-Znajek mechanism is working. However, we conclude that the jet is launched mainly by 
the magnetic field amplified by the gravitational collapse and differential rotation around the black 
hole, rather than the Blandford-Znajek mechanism. 

Subject headings: gamma rays: bursts — relativity — black hole physics — accretion, accretion discs 
— supernovae: general 



1. INTRODUCTION 

Gamma-Ray Bursts (GRBs; in this study, we consider 
only long GRBs, so we refer to long GRBs as GRBs 
hereafter) have be en mysterious phenom ena since their 
discovery m 1969 (jKlebesadel et al. 19731 ). Last decade, 
observational evidence for supernovae (SNe) and GRBs 
association has been reported (e.g. Woosley and Bloom 
2006, and references therein). 

Some of the SNe that associate with GRBs were very 
energetic and blight. The estimated explosion energy was 
of the order of 10 52 ergs, and produced nickel mass was 
~ 0.5M©. Thus they are categorized as a new type of 
SNe (sometimes called as hypernovae). The largeness of 
the explosion energy is very important, because it can not 
be explained by the standard core-collapse SN scenario, 
and other mechanism should be working at the center of 
the progenitors. 

The promising sce narios are the collapsar sce- 
nario dWooslev 19931) and the magnetar sce- 
nario (jUsov 19921 ). In the collapsar scenario, a 
rapidly rotating black hole (BH) is formed at the 
center, while a rapidly rotating neutron star with strong 
magnetic fields (~ 10 15 G) is formed in the magnetar 
scenario. Many numerica l simulations have been done 
for the collapsar s cenario dMacFadv en fc Wooslev 19991 : 
Proga et al. 20031: iProp fc Beeelman 2003] 



and the magnetar scen ario dTakiwaki et al. 2004 
Komissarov fc Barkov 20071 \B urrows et al. 2 007: 

Bucciantini et al. 20081 : Dessart et al. 2008 ; 

Takiwaki et al. 20081 iBucciantini et al. 20091 ) ~ In 
this study, we investigate the collapsar scenario. 

In the collapsar scenario, a BH is formed as a result 
of gravitational collapse. Also, rotation of the pro- 
genitor plays an essential role. Due to the rotation, 
an accretion disk is formed around the equatorial 
plane. On the other hand, the matter around the 
rotation axis falls into the BH almost freely. It is 
pointed out that the jet-induced explosion along the 
rotation axis may occur due to the heating through 
pair annihilation of neutrinos and anti- neutrinos that 
are emitted from the accr et ion disk dWooslev 19931 : 
iMacFadven fc Wooslev 19991 : iFrver fc Meszaros 20001) . 
Effect of extraction of rotation energy from the accretion 
disk by magnetic field lin es that leave the disk sur face 
(Blandford-Payne effect (jBlandford fc Pavne 1982])) is 
also investigated by seve ral autho rs dProga et al. 2003 : 



Proga fc Beeelman 20031: 



IMizuno et al. 2004a 



Mizuno et al 2004at IMizuno et al. 2004b| 

Proea 20051: IFuiimoto et al. 2006t iShibata et al. 2006 1 
Nagataki et al. 20071: ~ iSekiguchi fc Shibata~2 007 
Suwa et al. 2007t iBarkov fc Komissarov 2008aT) 
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Mizuno et al. 2004b!: iProga 20051: IFuiimoto et al. 2006 : 
Nagataki et al. 2007t iSuwa et al. 2007m . Recently, 
the effect of extraction of rotation energy from the 
BH through outgoing poynting flux (Blandford- 
Znajek effect (iBlandford fc Znaiek 19771) 1 is inves- 
tigated (|Barkov fc Komissarov 2008af) . In order to 
investigate the collapsar scenario completely, a high- 
quality numerical code including effects of a lot of 
microphysics (neutrino physics, nuclear physics, and 
equation of state for dense matter) and macrophysics 
(magneto-hydrodynamics, general relativity) has to be 
developed. Although many numerical studies have been 
reported, such a numerical code has not been developed 
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yet. Thus we have to develop our numerical code step 
by step. 

In this study, we investigate the dynamics of collapsars 
taking into account the general relativistic effects. Ex- 
traction of rotation energy from a rotating BH is one of 
them. Also, even when the rotation energy is extracted 
from the accretion disk, the properties of the accretion 
disk should depend on the properties of the BH: if the 
BH is rotating, the inner region of the accretion disk 
should be enforced to co-rotates with the BH. We inves- 
tigate how a jet is launched at the center of a progenitor, 
and how the property of the jet is. Effects of rotation 
of the BH on the formation of GRB jet have not been 
investigated so much. Barkov and Komissarov (2008) 
is a pioneering study. However, only one case is inves- 
tigated in their study, and the initial progenitor model 
they used is a simplified one- dimensional m odel without 
rotation and magnetic fields (Bcthe 1990). Since there 
should be many initial conditions of progenitors (pro- 
genitor mass, metallicity, angular momentum, magnetic 
fields), it should be important to investigate the general 
relativistic effects using a different initial condition from 
the previous study. In this study, we use a realistic initial 
condition for the progenitor model that is developed by 
Woosley and Heger (2006), in which rotation and mag- 
netic fields are taken into account. 

When we investigate the general relativistic effects, 
one has to develop a General Relativistic Magneto- 
Hydro Dynamic (GRMHD) code. So far, there are 
many studies on GRMHD code for fixed background 
space times using high-order conservative schemes based 
on either appr oximate or full wave- d ecomposition Rie- 
mann solvers (fGammie et al. 20031: iKomissarov 2005 : 
Anninos et al. 20051: lAnton et al. 2006 : 

Del Zanna et al. 2007t pTchekhovs kov et al. 20071) or 
non-conservative sc hemes (|De Villiers fc Hawlev 20031 : 
lAnninos et al. 20051 ). Since the accreted mass onto the 
BH is still less than the initial BH mass in this study, 
we take the GRMHD code for the fixed background. 
Especially, we develop our code using the conservative 
scheme of Gammie et al. (2003) with the method of 
Noble et al. (2006) for transforming conserved variables 
to primitive variables. 

The plan of the paper is as follows. In section O we 
present the formulation of the GRMHD code. In sec- 
tion^ we show results of many, well-known test calcula- 
tions to confirm the reliability of the code. After we show 
the reliability, we present results of numerical simulations 
of collapsars in section 0J Summary and discussion are 
presented in section 

2. DEVELOPMENT OF GRMHD CODE 

We have developed a two-dimensional GRMHD code 
following Gammie et al. (2003) and Noble et al. 
(2006). We have adopted a conservative, shock-capturing 
schem e with Harten, La x, and van Leer (HLL) flux 
term (H arten et al. 1983T ) wit h flux-inte rpolated con- 
strained transport technique (|T6th 2000D . We use a 
third-order Total Variation Diminishing (TVD) Runge- 
Kutta method for evolution in time, while monotonized 
central slope-limited linear interpola tion method is used 
for second-order accuracy in space (|van Leer 19771) . 2D 
scheme (2-dimensional Newton-Raphson method) is usu- 
ally adopted for transforming conserved variables to 



primitive variables (|Noble et al. 2006T ). 

When we perform simulations of GRMHD, Modified 
Kerr-Schild coordinate is basically adopted with mass 
of the BH (M) fixed where the Kerr-Schild radius r is 
replaced by the logarithmic radial coordinate x\ = In r. 
When we show the result, the coordinates are sometimes 
transfered from Modified Kerr-Schild coordinate to Kerr- 
Schild one for convenience. In the following, we use G = 
M = c = 1 unit. G is the gravitational constant, c is the 
speed of light, and M is the gravitational mass of the 
BH at the center. Throughout this pa per we follow the 
standard notation (M isner et al. 1970 ). 

2.1. Formalism 

Number of variables that appear in the equations of 
GRMHD is 13: rest-mass density (p), internal energy 
density (u), pressure (p), four- velocity of fluid (u M ), 
and Faraday tensor (F^ v ). Note that Faraday tensor 
has only 6 independent components due to the relation 
F^" = —F Ufi . We can reduce the number of independent 
variables to 8 using the MHD condition {u^F^" = 0), 
equation of state (p = (7 — l)u: 7-law gas is assumed), 
and the unit length of the four velocity (u^u^ = — 1). 
Note that the number of independent equations of the 
MHD condition is 3. We choose (/3,w,m%B 1 ) as the 8 in- 
dependent variables where u 1 is the space component of 
the four velocity. B % can be written as Sft 1 /a where a is 

the lapse function (a = an( i ^ 1S ^ ne magnetic 

field measured by the Fiducial observer (FIDO) whose 
four velocity is — (—a, 0,0,0). We call these inde- 
pendent variables as the primitive variables. Below, we 
introduce the conserved variables. Of course, number of 
the conserved variables is also 8. Thus we require 8 basic 
equations to follow the time evolution of the system. 

The basic equations of GRMHD represent the rest- 
mass conservation, the energy-momentum conservation, 
and space component of the induction equation that de- 
termines the time evolution of the magnetic fields. These 
are: 

dt^gpu 1 ) = -&(V=spu*) (l) 
dt{V~gTi) = -bW^sFi) + y/=sWF™ (2) 

d t {^—gB l ) = -d 3 [V=$(&V - Vu*)] , (3) 

where T^ v is the stress energy tensor that is composed of 
the sum of the matter part (T^ attCT — (p + u+p)u^u u + 
pg^) and electromagnetic part (T^ = F^ a F£ - 
g^F a pF al3 /A). The factor of y/A ^ is absorbed into th e 
definition of the Faraday tensor (jGammie et al. 2003). 
6 M is introduced so that Eq.Q looks simple, and 
it is defined as = e^ KX u ly F\ K , where e» VKX = 
(—l/y/—g) [pv\n\. [p,v\k] is the completely antisymmet- 
ric symbol. In the fluid-rest frame, becomes (0,B l ). 

In this study, we adopt the conservative scheme for 
integration of the GRMHD equations. In this case, the 
left terms of Eq.(p}-((3]) are considered to be fundamental 
variables and called as the conserved variables. The right 
terms of Eq.([T])-((3]) are flux terms with a source term (the 
second right term of Eq.Q). 

Since we have to estimate pressure of the fluid, we 
have to estimate the primitive variables from the con- 
served variables at each time step. The problem is that 
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the primitive variables can not be expressed analytically 
by the conserved variables. Thus we have to use the 
Newton-Raphson method to o btain the primitiv e vari- 
ables from the conserved ones (jNoble et al. 20061) . 

Basically, we adopt the 2D scheme introduced by No- 
ble et al. (2006) to calculate the primitive variables. 
However, it sometimes happens that the 2D scheme fails 
to converge well, and the primitive variables can not be 
obtained precisely. In such a case, we first adopt the 
lDw scheme introduced by Noble et al. (2006) and see 
whether the lDw scheme converges. If it converges well, 
we adopt the primitive variables obtained by the lD\y 
scheme for the next time step. Otherwise, we adopt the 
second choice explained in the following subsection. 

2.2. Supplemental Method to Calculate Primitive 
Variables 

Following Noble et al. (2006), we introduce convenient 
variables v 2 , W, Q M , and Q M . These variables, of course, 
depend on the primitive variables. The definition of the_se 
variables are: v 2 = W = ujT 2 , = aT 1 ^, and 
— j\Q X > where if is the fluid velocity relative to FIDO, 
uj = p + u + p, T = 1/Vl - v 2 , and j^x = g^x+n^nx- It 
is apparent that and Q M can be written analytically 
by the conserved variables. On the other hand, v 2 and 
W can not be expressed analytically by the conserved 
variables. Thus, we have to solve v 2 and W numerically 
in order to determine the proper, corresponding primitive 
variables. 

Here we show that an upper limit and a lower limit 
for W can be obtained before searching for a solution of 
W and v 2 numerically. Thanks to this fact, all we have 
to do is to seek the solution with the condition W m i n ^ 
W ^ W max . From Eq.(28) and Eq.(29) in Noble et al. 
(2006), W and v 2 satisfy the following equations: 



y oq29 



t/ oq28 

2 



Q 2 W 2 + (Qf,W) 2 (^ 2 + 2W) 



2W 2 



(5ft 2 + W) 2 W 2 
-W + p-iQun") 



(4) 



1. (5) 



From these equations and the relation ^ v 2 < 1, v 2 
and W satisfy the following relations: 

f(W) = W 4 + 2$l 2 W 3 + (5R 4 - Q 2 )W 2 - 2(Q^W) 2 W 
-^(Q^Wf^O (6) 

g (W) = W 3 + { l -^ 2 + (QX) - pW 2 - \{Q^f ^ 

(7) 

h(W) = W 3 + {ft 2 + {Q^)-p}W 2 - ^{Q^Wf ^ 0. 

(8) 

Since /(0) ^ 0, / (0) 5^ 0, and at least one of the solution 
for / (W) — is less than 0, there is only one positive 
solution W a that satisfies f(W a ) — 0. Thus, from Eq.©, 
W has to be greater than W a . 

We can understand the behavior of g(W) from its first 
derivative for W: 





-1 o 1 





g (W) = w 



3W + 2{^ 2 + {Q^)-p} 



(9) 



Fig. 1. — Simulation of ID shock tube test (Komissarov 1999). 
The state at t = 1.0 is shown in the figure. Number of grid points 
is 600. The calculation region is set to be —2 < x < 2. The upper 
left panel shows density, the upper right panel shows pressure, the 
lower left panel shows the velocity in the x-direction, and the lower 
right panel shows the bulk Lorentz factor. 

It is apparent that W — is a solution for g (W) = 0. 
As for the other solution(s), it is not so obvious because 
the pressure p depends on W and v 2 . However, it will 
be natural to consider that the monotonic relation holds 
between W and p. It means that the pressure rises when 
W becomes larger. If this assumption is adopted, as 
long as g (W) = has another solution, it is a positive 
one W = W a = 0. This is because when W — 0, p 
should be also and [3W + 2{l/23? 2 + (Q^) - p}} is 
a positive value. Thus, g(W) = has only one positive 

solution Wb- This holds even if g (W) = has only one 
solution at W = 0. Also, same conclusion can be derived 
for h(W): there is only one positive solution W c that 
satisfies h(W c ) = 0. 

Since h(W) ^ g(W), the relation W c ^ W b holds. 
In conclusion, W has to be in the range W m [ n = 
Max(Vy a , W c ) S W ^ W b = W max . Thus all we have to 
do is to find a solution of W that satisfies v 2 q28 = v 2 q29 
in this range. This procedure is more expensive than the 
2D scheme and the lDw scheme, but the solution for W 
and v 2 is more likely to be found because the range for 
the solution of W is determined apriori. Thus we use this 
method as a supplementary one to obtain the primitive 
variables. 

3. TEST CALCULATIONS 

Using the GRMHD code that is developed in this 
study, we check whether it can pass many, well-known 
test calculations. The first three tests are special rel- 
ativistic hydrodynamic (SRHD) or special relativistic 
magnetohydrodynamic (SRMHD) calculations, while the 
rest of three tests are GRMHD ones. 

3.1. Shock Tube Problems 



4 



S. Nagataki 



30 
20 
10 




M 





I I I I I 

: n 


i i 

By : 






', , , , i , , , 





-1 1 



Fig. 2. — Simulation of ID collision test (Komissarov 1999). The 
state at t = 1.2 is shown in the figure. Number of grid points is 600. 
The calculation region is set to be —2 < x < 2. The left panels 
show density, velocity in the x-direction, and bulk Lorentz factor 
(from top to bottom), while right panels show pressure, velocity 
in x-direction, and y-component of magnetic field (from top to 
bottom). 

ID shock tube tests are the most basic test problems 
for SRHD/SRMHD. We have carried out a number of 
the test simulations introduced in Komissarov (1999) 
and Balsara (2001). Here we describe only two of them. 
One is the shock tube testl and the other is the collision 
test (jKomissarov 19991: iMizuno et al. 20061) . The initial 
left and right states are summarized in Tabled] Number 
of grid points is 600 for both simulations. The results 
are shown in FigQ] and Fig[2l which show that the test 
calculations are well solved as in the previous studies. 

3.2. Double Shock Problems 

Here 2D shock tube problem is done to con- 
firm whether the shock dynamics in the multidimen- 
sional flow can be solved safely. This problem in- 
cludes the interactions of shocks, rarefactions, con- 
tact discontinuities. Initially a square computational 
domain is prepared in x-y plane and divided into 
four quarter boxes. Initial condition in each box 
is summarized in Table Ql This condition is same 
with previous study (iDel Zanna fc Bucciantini 20011 : 
iZhang fe MacFadven 20061 : iMizuta et al. 20061 ). We use 
400x400 uniform grid points in a square computational 
box. Boundary condition is open ones. Density contour 
at the final stage of the simulation is shown in Fig(3j 
which shows that our code reproduces the previous stud- 
ies very well. 

3.3. Cylindrical Explosion Test 

Here we go to a SRMHD tes t. A famous, cylindri- 
cal blast explosion te st is done (|Del Zanna et al. 20031 : 
iLeismann et al. 20051) . We use the [0,1] x [0,1] Carte- 




sian grid with a resolution of N x 



250 grid 



Fig. 3. — Simulation of 2D shock tube problem. Density contour 
at t = 0.4 is shown in the figure. Numbers of grid points are 
400x400. The calculation region is set to be < x < 1 and 
< y < 1. 

points. We define an initially static background with 
p = 1.0, p = 0.01, and B x = 4.0. The relativistic flow 
comes out by setting a much higher pressure, p = 10 
within a circle of radius r = 0.08 placed at the center of 
the domain. 7 for the equation of state is set to be 4/3. 
Final time is set to be 0.4. The result is shown in FigJH 
The upper left panel shows the density contour in loga- 
rithmic scale. The upper right panel shows the pressure 
contour in logarithmic scale. The lower left panel shows 
contour of the bulk Lorentz factor. The lower right panel 
shows the divergence of the magnetic fields in logarithmic 
scale with magnetic field lines. These results are consis- 
tent with the previous studies. Especially, the divergence 
of the magnetic fields is kept as small as 10~ 14 . 

3.4. Gammie's Flow 

Next we consider a GRMHD test. A steady, 
magnetized inflow solution on the equatorial plane 
around a Kerr BH is considered (|Takahashi et al. 19901 : 
IGamm ic 1999). Initially, the steady inflow solution for 
the Kerr parameter a — 0.5 and the magnetization pa- 
rameter Fg^ — 0.5 is set, and time evolution of the sys- 
tem is followed by the GRMHD code. In this calculation, 
Boyer-Lindquist coordinate is used. The calculation re- 
gion is set to be [2.0 ^ r ^ 4.04] and [0.5 - 10~ 3 ^ 
0/tt ^ 0.5 + 10" 3 ]. The model is run for t = 1.5. The 
physical values at boundaries are fixed throughout the 
simulation. Results are shown in FigO density, ra- 
dial component of the 4-velocity, the <j) component of 
the 4-velocity, and at the final stage of the simu- 
lation. When the initial state is written in the same 
figure, we can see that the final state coincides with 
the initial state. To show it more quantitatively, we 
introduce the norms of the errors for these values as 
a function of the number (N) of grid points in the ra- 
dial coordinate. The definition of the norm of the error 
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Fig. 4. — A RMHD 2D cylindrical explosion test with a pressure 
jump as high as 10 5 . The state at t = 0.4 is shown in the figure. 
Numbers of grid points are 250x250. The calculation region is set 
to be < x < 1 and < y < 1. The upper left panel shows the 
density contour in logarithmic scale. The upper right panel shows 
the pressure contour in logarithmic scale. The lower left panel 
shows contour of the bulk Lorentz factor. The lower right panel 
shows divergence of the magnetic fields in logarithmic scale with 
magnetic field lines. 

is £*zf |a(fmal) - a(initial)| /Efzf | a (initial) j . In FigH 
the norms of errors are shown. We can see that these val- 
ues converge roughly proportional to N~ 2 , as expected. 

3.5. Blandford-Znajek Monopole Solution 

Further we continue to test the GRMHD code. 
We consider the Blandford-Z najek monopole solu- 
tion (Blandfor d fc Znai ek 1977). This analytic solu- 
tion ha s been investigated numerically by previous 
studies (iKomissarov 2004b: McKinnev & Gammie 2004; 
iTanabe fc Nagataki 2008f T 

The computational domain is axisymmetric, with a 
grid that extends from r- m = 0.98r + to r out = 230 and 
from 9 = to = n where r + = 1 + \J\ — a 2 is the outer 
event horizon. The numerical resolution is 300 x 300. As 
an initial condition, we put the 0t h order terms of th e 
monopole solution around the BH (Komissarov 2004b). 
That is, W = -nlF^ v = (0, asin^/y^, 0, 0) in the 
Kerr-Schild coordinate where *F^ V and g are the dual 
field tensor and determinant of the Kerr-Schild met- 
ric. The plasma velocity relative to the FIDO is set 
to zero initially, and its pressure and density are set 
to small value (P = p = 5R 2 /100) so that the sys- 
tem becomes Force-Free like. Also, to keep the mag- 
netization reasonably low, when the critical condition 
0MB 2 > T 2 p + ( 7 r 2 - (7 - l))u is satisfied, density 
and internal energy are increased by the same facto r so 
that the critical condition holds (jKomissarov 20040 ). 7 
is set to be 4/3. We have performed numerical simula- 
tions with the Kerr parameters 0, 0.01, 0.1, 0.2, 0.3, 0.4, 
0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, and 0.995 until time t = 
200. 

The total energy flux, which is the integrated outgoing 
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Fig. 5. — Gammie's equatorial inflow solution in the Kerr metric 
with a = 0.5 and magnetization parameter Fg^ = 0.5. Number 
of grid point is 1024. The state at t = 1.5 is shown in the figure. 
The panels show density, radial component of the 4-velocity, the 
4> component of the 4-velocity, and SR^ at the final stage of the 
simulation. Boyer-Lindquist coordinate is used for the simulation. 
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10 100 1000 
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Fig. 6. — Convergence results for the Gammie's equatorial in- 
flow solution in the Kerr metric with a = 0.5 and magnetization 
parameter Fg^ = 0.5. Norms of the error for p, u r , , and 
at the final stage of the simulation are shown in the figure. The 
straight line represents the slope expected for second-order conver- 
gence. The definition of the norm of the error is written in the 
text. 
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Fig. 7. — Outgoing poynting fluxes as a function of the zenith an- 
gle for the Blandford-Znajek monopole solution with Kerr param- 
eter a = 0.01,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, and 0.95. 
The fluxes are measured at r = 20 and t = 200. Numbers of grid 
point are 300x300. 

poynting flux over the zenith angle, can be written as 



E = 2tt 



-g(-T[) = 2n dx 2 — 



2tt / dx 2 F E , 




wher e x 2 = 9/tt is intr oduced as a convenient vari- 
able (jGammie et al. 20031 ). 

In Fig|7l the outgoing poynting fluxes (Fe) as a func- 
tion of zenith angle are shown. The fluxes are measured 
at r = 20 at the final stage of the simulations. We would 
like to note that the outgoing poynting flux hardly de- 
pends on the radius where it is evaluated. This means 
that the conservation of the outgoing poynting flux is 
confirmed numerically. 

In FigJSJa), we plot the total energy flux (E) at the fi- 
nal stage for small Kerr parameters (0 < a < 0.2) by rect- 
angular points. Dashed line is just the interpolation of 
the calculated values. For comparison, the second-order 
analytical solution is shown by dotted line and the forth- 
order analytical solution is shown by solid line. From 
this comparison, we can see that all of them coincide 
with each other. Thus the results of the numerical simu- 
lations by the GRMHD code are confirmed by analytical 
solutions. 

The situation becomes different for large Kerr parame- 
ters. In FiglUJb), we plot the same values with Fig. EJa), 
but for wide range of the Kerr parameters (0 < a < I). 
We can see clearly the difference among three cases. This 
is because the analytical solution is obtained by the per- 
turbation method in Kerr parameter, and it is applicable 
only for small Kerr parameters. Of course, there is no 
such limitation for the numerical simulations. Thus the 
total energy flux obtained by the numerical simulation is 
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Fig. 8. — Upper panel (a): Comparison of the derived, conserved, 
total energy flux. Dashed line with rectangular points is numerical 
result for small Kerr parameter (0 < a < 0.2), dotted line shows 
the second-order analytical solution, and solid line represents the 
forth-order analytical solution. Lower panel (b): Same with upper 
panel, but for wide range of the Kerr parameters (0 < a < 1). 
Simulations are done for the Kerr parameters 0, 0.01, 0.1, 0.2, 0.3, 
0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, and 0.995 until t = 200. 

more reliable than the analytical estimation (see Tanabe 
and Nagataki (2008) for detailed discussion). 

3.6. Fishbone and Moncrief's Test 

Here we present a final te st of the GRMHD code. A 
steady and stationary to rus ([Fishbone Sz Moncrief 19761 : 
lAbramowicz et al. 19781 ) around a Kerr BH that is sup- 
ported by both centrifugal force and pressure is solved 
numerically. Of course, it should be solved as a steady 
and stationary state. 

We have integrated a Fishbone-Moncrief solution 
around a Kerr BH with a = 0.9. We set w'u^ = 4.45 and 
B: m = 6.0. The grid extends radially from r- ul — 1.40 to 
Tout = 100. The same floors with Gammie et al. (2003) 
are used for p and u. The numerical resolution is N x N 
and the solution is integrated for t = 10. The resulting 
norm of the error, which converges roughly proportional 
to A^ 2 , is shown in FigJ5] 

Next we follow the time evolution of the Fishbone- 
Moncrief solution with magnetic fields. The vector po- 
tential, Afi cx max(p/p max — 0.2, 0) where p max is the pea k 
density in the torus, is introduced (jGammie et al. 2 003). 
The field is normalized so that the minimum value of 
Pgas/Pmag becomes 10 2 . The time integration extends for 
t = 2000. The number of grid points is 256 x 256, and 
the grid extends radially from r in = 1.40 to r out = 300 
while it extends in the zenith angle from 9 — to 8 = ir. 

The density contours in logarithmic scale (from 10 -6 
to 10 3 ) are shown in Fig|TU] These are projected on the 
(r sin f?,r cos f?)-plane. The upper left panel shows the 
initial state. The upper right panel shows the final state 
of the simulation without magnetic fields. The lower left 
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Fig. 9. — Convergence results for the Fishbone-Moncrief problem 
for a Kerr BH with a = 0.9. The straight line represents the slope 
expected for second-order convergence. Kerr-Schild coordinate is 
used for the simulation. 

panel shows the final state of the simulation with mag- 
netic fields. The lower right panel is same with the lower 
left one, but for a wide region. Due to the presence of 
the magnetic fields, the angular momentum in the torus 
is conveyed outward and the the torus starts to accrete, 
and the jet is launched from the BH around the polar 
regio n. This result is con sis tent with the previous stud- 
ies (iGammie et al. 20031 IMcKinnev k. Gammie 2004; 
IMcKinnev 2006al ; iMckinnev 2006H ). 

4. SIMULATION OF COLLAPSAR 

Since our code has passed the many test calculations 
shown in the previous section, we now simulate the dy- 
namics of a collapsar using the code. However, we have 
to say beforehand that no microphysics is included in 
the code such as nuclear reactions, neutrino processes, 
and equation of state for dense matter. So this is the 
FIRST STEP of our project to simulate the dynamics of 
a collapsar and formation of a relativistic jet of a GRB. 

4.1. Method of Calculation 

We have done a 2D GRMHD simulation of a col- 
lapsar using the Modified Kerr-Schild coordinate and 
G = c = M = 1 units. When we show results, the 
coordinate is transfered from the Modified Kerr-Schild 
coordinate to the Kerr-Schild one, and the units are fre- 
quently transfered to cgs units. The calculated region 
corresponds to a quarter of the meridian plane under the 
assumption of axisymmetry and equatorial symmetry. 
The spherical mesh with 256(r)x 128(6*) grid points is 
used for all the computations. The calculated region cov- 
ers from r =1.8to3xl0 4 (that corresponds to5.3xl0 5 cm 
and 8.9xl0 9 cm in cgs units, as explained below) with 
uniform grids in the Modified Kerr-Schild space. 

We adopt the model 12TJ in Woosley and Heger 
(2006). This model corresponds to a star that has 12M 




Fig. 10. — Density contour in logarithmic scale (from 10 6 to 
10 3 ) for the Fishbone-Moncrief problem for the Kerr parameter 
a = 0.9. Number of grid points is 256 X 256. The simulations 
are done until t = 2000. The upper left panel shows the initial 
state. The upper right panel shows the final state of the simulation 
without magnetic fields. The lower left panel shows the result with 
magnetic fields. The lower right panel is same with the lower left 
one, but for a wide region. These results are projected on the (r 
sin 9,i cos 6>)-plane. 

initially with 1% of solar metallicity, and rotates rapidly 
and does not lose its angular momentum so much by 
adopting small mass loss rate. As a result, this star 
has a relatively large iron core of 1.82M , and rotates 
rapidly (the estimated Kerr parameter that a BH forming 
of mass and angular momentum of the inner 3 M would 
formally have is 0.57) at the final stage. Of course, what 
kind of stars are ap propriate for prog enitors of GRBs is 
still under debate (jYoon et al. 2 006). Thus we chosed 
the model 12TJ as a first example of our study because 
the iron core is large and rotating rapidly, which seems 
to form a rapidly-rotating BH, among the models listed 
in Woosley and Heger (2006). We assume that the cen- 
tral part of the star with 2M has collapsed and formed 
a BH at the center with the Kerr parameter a = 0.5. 
We also assume that the gravitational mass of the BH is 
unchanged throughout the calculation. Since M = 2M , 
r = 1 corresponds to 2.95 xl0 5 cm, as explained above. 
Also, the inner boundary r = 1.8 is set within the outer 
horizon r + = 1 + \/l — a? = 1.866. 

Since 1-D calculation is done for the model 12TJ, we 
can use the data directly only for the physical quanta 
on the equatorial plane. As for the density, internal en- 
ergy density, and radial velocity, we assume the struc- 
ture of the star is spherically symmetric. We also set 
u e = initially. As for u^, we extrapolate its value such 
as u^(r, 9) = u*{r, tt/2) x sinfl. 

Effects of magnetic fields are taken into account in 
the model 12TJ. However, again, since 1-D calculation is 
done, we do not know the configuration of the magnetic 
fields. It is difficult to extrapolate magnetic fields that 
satisfy the condition divB = everywhere. Also, there 
are much uncertainty on the amplitude of the magnetic 
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Fig. 11. — Contours of rest mass density at the central region 
in logarithmic scale, in which cgs units are used assuming that 
the gravitational mass of the BH is 2Mq . The length unit in the 
vertical/horizontal axes corresponds to 2.95 xlO 5 cm. Upper panel 
(a) shows the state at t = 110000 (that corresponds to 1.0835 sec), 
while lower panel (b) shows the one at t = 180000 (that corresponds 
to 1.773 sec). These results are projected on the (r sin 0,r cos 6)- 
plane. 

fields in a progenitor. Thus we do not use the informa- 
tion on magnetic fields of the model 12TJ. Rather, we 
adopt the same treatment in section 13.61 That is, the 
vector potential oc max(p/p max — 0.2,0) sin 9 where 
Pmax is the peak density in the progenitor (after extract- 
ing the central part of the progenitor that has collapsed 
and formed a BH). The field is normalized so that the 
minimum value of p ga s/Pmag becomes 10 2 . The defini- 
tion of Pmag is Pmag = b 2 /2. The reason why we adopt 
the strong dependence on the zenith angle for A^, is so 
as not to suffer from discontinuity of magnetic fields at 
the polar axis. The resulting biggest amplitude of the 
magnetic fields is 7.4xl0 8 G at r = 950 (2.8xl0 8 cm). 

We use a simple equation of state p gas = (7 — V)u 
where we set 7=4/3 so that the equation of state roughly 
represents radiation gas. 

As for the boundary condition in the radial direction, 




Fig. 12. — Same with Fig. 11(b), but for a wider region. 



we adopt the outflow boundary condition for the inner 
and outer boundaries (|Gammie et al. 20031 ). As for the 
boundary condition in the zenith angle direction, axis 
of symmetry condition is adopted for the rotation axis, 
while the reflecting boundary condition is adopted for the 
equatorial plane. As for the magnetic fields, the equato- 
rial symmetry boundary condition, in which the normal 
component is continuous and the tangential component 
is reflected, is adopted. 

4.2. Results 

In FigfTTJ color contours of rest mass density at the 
central region are shown. Colors represent the density 
in units of g cm~ 3 in logarithmic scale. These results 
are projected on the (r sinf?,r cos 6*)-plane. The length 
r = 200 corresponds to 5.9 xlO 7 cm. The time unit 
corresponds to 9.85 xlO -6 sec. Upper panel (a) repre- 
sents the contours of rest mass density at t = 110000 
(that corresponds to 1.0835 sec), while lower panel shows 
the contours at t = 180000 (that corresponds to 1.773 
sec). FigfTUis the same figure with FigfTlTb'). but for a 
wider region. A jet is clearly seen along the rotation axis. 
In Fig U31 mass accretion rate history on the horizon is 
shown. The definition of the mass accretion rate is 

M = 2 x 2tt \ dOy/^jpu 7 '. (11) 
Jo 

It takes about 0.15 sec for the inner edge of the matter 
to reach the horizon. When the matter reaches there, 
there is an initial spike of the mass accretion rate. After 
that, there is a quasi-steady state like Fig fTTT a) is real- 
ized. Then, the jet is launched at ~1.1 sec. After that, 
the mass accretion rate varies rapidly with time, which 
resembles to a typical time profile of a GRB. 

We show color contours of the plasma beta (pg as /p mag ) 
in logarithmic scale at t = 180000 in Figll4l As ex- 
pected, the plasma beta is low in the jet region while it 
is high in the accretion disk region. We show color con- 
tours of bulk Lorentz factor around the central region 
at t = 180000 in FigfTBTa) (upper panel, in logarithmic 
scale). Color contours of the energy flux per unit rest- 
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Fig. 13. — Mass accretion rate history on the horizon. The unit 
Mq sec -1 is used assuming that the gravitational mass of the BH 
is 2Mq throughout the calculation. 

mass flux (E = —T[/(pu r )), which is conserved for an 
inviscid fluid flow of magnetized plasma, are also shown 
in Fig[T5Tb') (lower panel, in logarithmic scale). This 
value represents the bulk Lorentz factor (Too) of the in- 
vischid fluid element when all of the internal and mag- 
netic ener gy are converted into kinetic energy at large 
distances (McKinncy 2006a). We can see that the bulk 
Lorentz factor of the jet is still low (Fig fTST a)). but it can 
be potentially as high as 10 2 at large radius (FigfloTb)). 
At t = 180000, the strength of the magnetic field (y/^b 2 ) 
at the bottom of the jet is found to be ~ 10 15 G, and 
u^/u* is ~ 0.1 at r ms on the equatorial plane. Here r ms 
is the marginally stable orbit. For the Kerr BH with 
a = 0.5, r ms is 4.23. As stated in section 14.11 the ini- 
tial biggest amplitude of the magnetic fields is 7 Ax 10 8 G 
at r — 950 where the initial density is ~ 10 6 g cm~ 3 , 
the expected amplification factor of the magnetic fields 
due to the gravitational collapse and differential rota- 
tion around the BH is (p/p ) 2/3 x (d^/dt/2-K) * At ~ 
100 x 0.016 x 180000 - 3 x 10 5 . Thus the initial mag- 
netic field can be amplified as large as several times of 
10 14 G, which is roughly consistent with the amplitude 
of the magnetic fields at the bottom of the jet. At late 
phase, the magneto-rotational instability (MRI) may be 
also working, which is discussed in the next section. 

In Fig |161 contours of the <f> component of the vec- 
tor potential {A^) at t = 180000 are shown. Level 
surfaces coincide with poloidal magnetic field lines, and 
field line density corresponds to poloidal field strength. 
As expected, the magnetic fields are strong at the jet 
region, which makes the plasma beta very low. From 
FigEl the opening angle of the jet is estimated as 
5° — 6°. From Fig fT4l [l"5T b) . andfTol this j et should corre - 
spond to the poynting flux jet (Hawley fc Krolik 20061 ). 
This jet is surrounded by the funnel- wall jet re- 




FlG. 14. — Contour of the plasma beta (p gas /pmag) at t = 180000 
in logarithmic scale. 



below. 

In Table [3] and [U the integrated energies of matter and 
electromagnetic field at t = 180000 are shown. The in- 
tegrated region is between the horizon and r = 200 (for 
Table[3l) or r = 40 (for Tabled, and within the zenith 
angle measured from the polar axis. As for the matter 
component, the contribution of the rest mass energy is 
subtracted. That is, 



E 



Matter 



= 2x2tt 



r=200 or 40 




Matter" 



-pu°u ). 
(12) 



Factor 2 is coming from the symmetry of the system with 
respect to the equatorial plane. The field part can be 
written as 



E 



EM 



4-7T 



r=200 or 40 



dr 



"S^ChEM- 



(13) 



It can be seen that the energy in electromagnetic field 
dominates that in matter within r < 40, while they be- 
come comparable within r < 200. Also, the integrated 
energy is still less than the typical explosion energy of a 
GRB (~ 5 x 10 50 erg; Frail et al. 2001). 

Finally, we show the rest-mass density, outgoing mass 
flux, and outgoing poynting flux in Fig ll7l The top 
panel (a) shows the rest-mass density (g cm~ 3 ) as a 
function of the zenith angle at r = 10r ms = 42.3 and 
t = 180000. It is seen that the low-density region 
is realized around the polar axis, which corresponds 
to the jet region (0.1 radian corresponds to 5.7°). 
The middle panel (b) shows the outgoing mass flux 
pu r (g cm" 2 s" 1 ) at r = 10r ms and t = 180000. It 
is seen that the outgoing mass flux exists around 
0.06 < 9 < 0.14, which corresp o nds to the funnel- 
wall jet dDe Villiers et al. 2003t iHirose et al. 2004 



McKinncy & Gammie 2004; iKato et al. 2004 

De Villiers et al. 20051: lHawlev fc Krolik 20061 ). The 
bottom panel (c) shows the outgoing poynting flux (in 
units of 10 50 erg s _1 rad" 1 ) at r = r + and t = 180000. 
Definition of the outgoing poynting flux is 



gion (|Hawlev fc Krolik 2006f ). which is shown in FigfTTl 



Fy 
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"S^r.EM- 
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Fig. 15. — Upper panel (a): contours of bulk Lorentz factor 
around the central region at t = 180000. Lower panel (b): con- 
tours of the energy flux per unit rest-mass flux at t = 180000 that 
represent the bulk Lorentz factor (Too) of the invischid fluid ele- 
ment when all of the internal and magnetic energy are converted 
into kinetic energy at large distances. The contours are written in 
logarithmic scale. 

Factor 2 is coming from the assumption of the symmetry 
of the system with respect to the equatorial plane. We 
can see that the positive outgoing poynting flux exists at 
the jet region (0 < 8 < 0.23 in radian). The integrated 
energy of the outgoing poynting flux at r = r+ and t ~ 
180000 is 4.6 x 10 46 erg s _1 . Since the duration of the 
jet in this study is ~ 0.7 s, this outgoing poynting flux 
seems to be too weak to explain the energy of the jet 
listed in Table [3] and QJ Thus we conclude that the jet is 
launched mainly by the magnetic field amplified by the 
gravitational collapse and differential rotation around the 
BH, rather than the Blandford-Znajek mechanism in this 
study. 

5. SUMMARY AND DISCUSSION 

In order to investigate the formation of relativistic jets 
at the center of a progenitor of a GRB, we have devel- 
oped a two-dimensional GRMHD code. In order to con- 
firm the reliability of the code, we have shown that the 
code passes many, well-known test calculations. Then 
we have performed a numerical simulation of a collap- 
sar using a realistic progenitor model. We have followed 




Fig. 16. — Contours of the component of the vector potential 
(A,p) at t = 180000. Level surfaces coincide with poloidal mag- 
netic field lines, and field line density corresponds to poloidal field 
strength. The biggest amplitude of the magnetic fields is ~ 10 15 G. 

the time evolution of the system for 1.773 s, and it was 
shown that a jet is launched from the center of the pro- 
genitor. We also found that the mass accretion rate is 
in quasi-stable state before the launch of the jet, while it 
shows rapid time variability that resembles to a typical 
time profile of a GRB after the launch. The structure 
of the jet is similar to the previous study: a poynting 
flux jet is surrounded by a funnel- wall jet. Even at the 
final stage of the simulation, the bulk Lorentz factor of 
the jet is still low, and the total energy of the jet is still 
as small as 10 48 erg. However, we found that the en- 
ergy flux per unit rest-mass flux [E = —T[ / (pu r )) is as 
high as 10 2 at the at the bottom of the jet. Thus we 
conclude that the bulk Lorentz factor of the jet can be 
potentially high when it propagates outward. Also, as 
long as the duration of the activity of the central engine 
is long enough, the total energy of the jet can be large 
enough to explain the typical explosion energy of a GRB 
(~ 5 x 10 50 erg). It is shown that the outgoing poynting 
flux exists at the horizon around the polar region, which 
proves that the Blandford-Znajek mechanism is working. 
However, we conclude that the jet is launched mainly by 
the magnetic field amplified by the gravitational collapse 
and differential rotation around the BH, rather than the 
Blandford-Znajek mechanism in this study. 

Whe n we apply the Bland ford-Znajek for- 
mula (jBarkov fc Komissarov 200 8b). the integrated 
outgoing poynting flux is 

£ = 3.6 x 10 50 /(a)*27 M 2~ 2 er S s_1 ( 15 ) 

where M 2 = M BH /2M = 1, * 2 7 = */10 27 G cm 2 , 
and /(a) = a 2 /(l + \/l - a 2 ) 2 = 0.07179. This value 
becomes 2.3 x 10 46 _B 2 5 erg s _1 for the jet with opening 
angle 9 = 5°, which is comparable to our numerical result 
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Fig. 17. — The top panel (a): the rest-mass density (g cm -3 ) as a 
function of the zenith angle at r = 10r ms = 42.3. The middle panel 
(b): the outgoing mass flux pu r (g cm -2 s — 1 ) at r = 10r ms . The 
bottom panel (c): the outgoing poynting flux (in units of 10 50 erg 
s -1 rad" 1 ) at r = r+. These values are evaluated at the final stage 
of the simulation (t = 180000). 



(Figure)). 

As for the efficiency of converting the released gravi- 
tational energy to the jet's energy, it can be estimated 
as follows: the mass accretion rate is ~ 0.1M© s _1 
(Fig[l3|), the total energy of the jet at the final stage 
is ~ 10 48 erg (TableOJ), and the duration of the jet is 
~ 0.7 s (Fig [T5|) . Thus the efficiency can be estimated 
as ~ 10 -5 . When we use the outgoing poynting flux at 
the horizon (4.6xl0 46 erg s -1 ), the efficiency is as low as 
~ 3 x 10~ 7 . These values seem to be very small com- 
pared with the previous stu dy (|De Villiers et al. 20051 : 
Mcki nnev fe Narav an 2007a]). One of the reason will be 
because they used an almost steady disk model. On the 
other hand, we used a realistic progenitor model that 
collapses gravitationally. Thus the resulting mass accre- 
tion rate is pretty high. Second reason may be because 
the efficiency is still low even at the final stage of the 
simulation. If we perform the simulation further, the ef- 
ficiency may become higher with time: mass accretion 
rate will become smaller, and the jet energy might be 
larger due to the amplification of the magnetic fields due 
to winding-up (and MM) effects. Also, when the initial 
amplitude of the magnetic field is set to be larger (as in 
Barkov and Komissarov 2008), the efficiency may be en- 
hanced. Further, we should investigate the dependence 
of the dynamics on progenitor models as well as the Kerr 
parameter of the BH. We are planning to investigate this 
point systematically in the next paper. 

It is well known that the system is unsta- 
ble against MRI when ther e is a strong nega- 
tive shear profile (dfl /dlnr) (Ba lbus fe Hawlev 19911 : 
iBalbus fe Hawlev"T9 94) , where f2 is the angular velocity. 
The saturation toroidal magnetic field strength is roughly 



expected to be B& ~ {Airp Y^rVL ()Akivama et al. 20031 : 
lAkivama fe Wheeler 2005|). which is confirme d by semi- 
global simulations (|Obergaulinger et al. 2 008) . The sat- 
uration poloidal magnetic field strength is roughly a n 
order of magnitude smaller (jObergaulinger et al. 2 006) . 
Thus may be amplified by MRI as strong as 1.4 x 
10 15 G pgr ms fl4 where is fi/10 4 rad s _1 . The char- 
acteristic timescale for saturating the MRI is the Alfven 

crossing time: tj± ~ 0.1ms Rep 9 B 15 where pg, 
and _Bi5 are the radius in units of 10 6 cm, the density in 
units of 10 9 g cm -3 , and the amplitude of magnetic fields 
in units of 10 15 G, respectively. Thus this characteristic 
timescale can be shorter than the winding-up timescale 
for strong magnetic fields. However, the length scale of 
the mode with the largest MRI growth rate is approx- 



-1/2 



where Pq.5 is the 



imately A M ri ~ 50 cm P . 5 B w p 9 
rotation period in units of 0.5 ms, which is too short 
to be resolved numerically. At least, it is not resolved 
in the beginning of the simulations. After the magnetic 
field is amplified to a certain value due to gravitational 
collapse an d winding-up effect, MRI may be working in 
this study (lOb crgaulin ger et al. 2006bl: fott et al. 20061 : 
iBurrows et al. 20071 : iDessart et al. 20081 ). It will be nec- 
essary to develop a sophisticated code that takes into 
account t he MRI effectively with hel p of semi-global sim- 
ulations (jObergaulinger et al. 20081 ) in order to evaluate 
the influence of MRI on the dynamics of a collapsar. 

It is well known that it becomes difficult to 
obtain the matter part of the primitive vari- 
ables (p ,u,u l ) precisely by the Newton-Raphson 
method (No ble et al. 20061) due t o numerical trun- 
cation errors (jKomissarov 200*21 lKomissarov"2 004a: 
Komissarov 2004bl iMcKinnev fc Gammie 2004 



Komissarov 20051: iMckinnev 20 06b)' when the elec 



tromagnetic part of the stress energy tensor ( T^) 
greatly exceeds the matter part ( ?M^ tter ). The problem 
is that the time integration of the electromagnetic part 
does not become so reliable, either. This is because the 
MHD condition (u^F^ = 0) is implicitly assumed in 
the basic equations, and the resulting basic equation 
of electromagnetic part depends on the velocity of 
fluid (Eqj3j. Such pathological conditions may be 
realized at the bottom of the jet in our study. In 
order to confirm the validity of our results in this 
study, we are planning to develop a general relativistic 
force-free code that is coupled w ith the GRMHD code 
sophisticatedly ([Mckinnev 2006^). 

It is very important to evaluate the terminal bulk 
Lorentz factor, because GRBs are considered to be 
emissions from relativistic flows with their bulk Loren tz 
factors greater than 10 2 (e.g. (jLithwick fc Sari 2001h ). 
Although an ad-hoc thermal (and kinetic) energy depo- 
sition into the polar region seems to l ead to relativistic 
jets with bulk L o rentz factors ~ 100 (lAlov et al. 2000 
Alov et al. 2002t iZhang et al. 20031: iZhang et al. 2004 



Cannizzo et al. 200 
Mizuta fc Alov 20081 : 
Wang et al. 2008h . it 



Miz uta et al. 2006 



Morsony et al. 2007 



such ad-hoc energy 
numerical simulations 
physics (Na gataki et al 
hand, numerical study on 
tromagnetically 



is still controversial whether 
deposition is justified by 
wit h proper neutrino 
20071) . On the other 

the acceleration of elec- 
powered jet requires quite high- 
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resolution (iKomissarov et al. 2007t iNaravan et al. 20071 : 
iTchekhovskov et al. 20"08[ IKomissarov et al. 20081) . 
Due to the reason, a simplified jet model with 
an idealized boundary condition is used at 
present in order to investigate whether the ini- 
tial poynting flux can be effectively converted 
into kinetic energ y (fTchckhovs kov et al. 20081 
IKomissarov et al. 2 008). According to their results, as 
long as confinement of the jet is realized, acceleration 
operates over several decades in radius and considerable 
fraction of the poynting flux can be converted into 
the kinetic energy. Thus, from the high-ratio of the 
poynting flux relative to rest-mass flux seen in our study 
(Fig[T5Tb)), a relativistic jet with high bulk Lorentz 
factor may be realized at large radius. 

It is true that the two-dimensional restriction can 
be a significant l imitation. First, anti-dynamo theo- 
rem (Moffat 1978) prevents the indefinite maintenance 
of the poloidal magnetic field in the face of dissipation. 
Second, axisymmetric simula tions tend to overem - 
phasize the channel mode (Haw lev fe Balbus 1992f ). 
which produces coherent internal magnetized flows 
rather than the more generic MHD turbulence. Hy- 
drodynamic instability in t he azimuthal direction ma; 
be also very important (|Nagakura fe Yamada 200 
iNagakura fc Yamada 20091 ). Thus we are plan- 
ning t o develop a three -dime nsional GRMHD 
code (iDe Villiers et al. 20031 : iHirose et al. 2004 : 
De Villiers et al. 20051 : lHawlev fc Krolik 2006 : 

Beckwith et al. 20081 : iShafee et al. 200S : 

and investigate the 
two-dimensional simulations of 



Mckinney & Blandford 2008) 



difference between 
collapsars with three-dimensional ones. 

In this study, we assumed that the central region 
of the progenitor has collapsed and a BH is formed 
at the center with surrounding envelope unchanged. 
Thus we solved the GRMHD equation on a fixed 
background. But the final goal of our project is to 
study how a GRB is formed from the gravitational 
collapse of a massive star. Thus we are planning to 
develop a GRMHD code on a dynamical background, 
which makes the study on the gravitational collapse 
and BH f ormation at t he center of a massive star 
possible (IShibata 2001 iSekigudii fc Shibata 2004 : 
Sekiguchi fc Shibata 20051: iBaiotti et al. 2005 : 

Duez et al. 2006HSekiguchi fc Shibata 20071 ). 

In this study, photo-disintegration of nuclei and 
neutrino processes are not taken into account. Photo- 
disintegration absorb considerable amount of thermal 
energy, and cooling/heating due to neutrino processes 
may have great influence on t he dynamics of a collap- 
sar dDi Matteo et al 20021: iKohri fc Mineshige 2002 1 
Nagataki et al 2003at iSurman fc McLaughlin 2004] 
Lee et al. 20051: IGu et al. 20061 : iNagataki et al. 200fl 
Kawanaka fc Mineshige 20071 : iKawabata et al. 2008 

iZhang fc Dai 2009 



Ros si et al. 20081 



Cannizzo fc Gehrels 20091) . 



Especially, 
neutrinos may 



pair- 
be a 



annihilation of electron-type 
key process to drive a GRB jet (TWooslev 1993 
MacFadycn & W oosley 19991 : E sano fc Fukuvama 20001 : 



Asano fc Fukuvama 200lt iMiller et al. 2003 : 

Surman fc McLaughlin 20051: |Kneller et al. 2006 : 

Shibata et al. 20071 : iBirkl et al. 2007f K Thus we are 
planning to include such microphysics in our code, and 
perform more realistic simulations of collapsars. 

The SNe associated with GRBs often show pe- 
culiar pr operties. Some are very energetic and 
blight (iGalama et al. 19981: llwamoto et al. 19981 : 
iHiorth et al. 20031 : iMalesani et al. 20041 ). but oth- 
ers prohi bit such blight SN e from being accom- 
panied dFvnbo et al. 20061 : iDella Valle et al. 20061 : 
iGal-Yam et al. 20061) . Since the brightness 

of SN e depends on the mass of produce d 
56 Ni (jWooslev et al. 19991 : iNakamura et al. 20011 ). 
it is suggested that there is a huge variety of the amount 
of 56 Ni in a SN that associates with a GRB. It is still 
controversial where and when 56 Ni is produced in a 
SN accompanied by a GRB (|Nagataki et al. 2003bl : 
INagataki et al. 20061). It m ay be produced in a 



GRB iet dMaeda et al. 20021: iMaeda fc Nomoto 2003 



Tanaka ct al. 2007; Macda et al. 2008; Tominaga 2009 



Maeda fc Tominaga 20091 : iBucciantini et al. 20091 ). or 



it may be produced in (o r outflow from) the accretion 
disk around th e BH dMacFadven fc Wooslev 1999 ; 
Pruet et al. 20041: iFuiimoto et al. 2004 : 

Surman et al. 20061 : IHu fc Peng 2008^ or it 
may be synthesized around a proto-neutron 
star (lUzdenskv fc MacFadven 20071 ). At present, it 
is impossible to investigate the explosive nucleosynthesis 
in a collapsar in our code because nuclear reactions are 
not taken into account. We are planning to include 
this effect, and study the site of 56 Ni production. 
Also, study of a GRB as a possible site where very 
heavy elements an d light eleme nt s are synthesized 
is very important dLemoine 20021 : IBeloborodov"200l 
iSuzuki fc Nagataki 20051 ). 
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S. Nagataki 



TABLE 1 

Initial Conditions for Shock Tube TestI and Collision Test 



Test Type 




P 


V 




V v 


V 




& 


B s 


Shock Tube Tcstl 


lelt state 


1 


1000 











1 










right state 


0.1 


1 











1 








Collision Test 


left state 


1 


1 


5/V26 








10 


10 







right state 


1 


1 


-5/V26 








10 


-10 






Note. — 7 for the equation of state is set to be 4/3. Final time is set to be 1.0 for Shock Tube Test 1, and 1.2 for Collision Test. 



TABLE 2 

Initial Conditions for 2D Shock Tube Problem 



Region 


X 


y 


P 


P 


v x 


yV 


A 


< x < 0.5 


0.5 < y < 1 


0.1 


1 


0.99 





B 


0.5 < x < 1 


0.5 < y < 1 


0.1 


0.01 








C 


< x < 0.5 


< y < 0.5 


0.5 


1 








D 


0.5 < x < 1 


< y < 0.5 


0.1 


1 





0.99 



Note. — 7 for the equation of state is set to be 5/3. Final time is set to be 0.4. 



TABLE 3 

Integrated Energies of Matter and Field (r<200) 






0.714" 


1.43 u 


2.14" 


2.86" 


3.57" 


Matter 
Field 


1.44E+46 
2.96E+46 


7.09E+46 
1.11E+47 


1.70E+47 
2.39E+47 


3.14E+47 
4.12E+47 


5.10E+47 
6.30E+47 


6 


4.29" 


5.00" 


5.71" 


6.42" 


7.14" 


Matter 
Field 


7.66E+47 
8.91E+47 


1.10E+48 
1.19E+48 


1.52E+48 
1.52E+48 


2.05E+48 
1.88E+48 


2.69E+48 
2.26E+48 



Note. — The energy is written in units of erg. As for the matter component, the contribution of the rest mass energy is subtracted. 



TABLE 4 

Integrated Energies of Matter and Field (r<40) 



e 


0.714" 


1.43° 


2.14" 


2.86" 


3.57" 


Matter 
Field 


6.89E+43 
8.60E+45 


3.15E+44 
3.44E+46 


8.96E+44 
7.73E+46 


2.03E+45 
1.37E+47 


3.95E+45 
2.14E+47 


e 


4.29° 


5.00° 


5.71" 


6.42" 


7.14" 


Matter 
Field 


6.76E+45 
3.08E+47 


1.04E+46 
4.19E+47 


1.47E+46 
5.46E+47 


1.96E+46 
6.91E+47 


2.50E+46 
8.52E+47 



Note. — Same with Table. 3, but integration is done for R<40. 
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